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Abstract 

Activation cascades are prevalent in cell signalling mechanisms. We study 
the classic model of linear activation cascades and find that in special but 
important cases the output of an entire cascade can be represented analyti- 
cally as a function of the input and a lower incomplete gamma function. We 
also show that if the inactivation rate of a single component is altered, the 
change induced at the output is independent of the position in the cascade 
of the modified component. We use our analytical results to show how one 
can reduce the number of equations and parameters in ODE models of cell 
signalling cascades, and how delay differential equation models can some- 
times be approximated through the use of simple expressions involving the 
incomplete gamma function. 

Keywords: Activation cascades, linear models, MAPK, gamma function, 
model reduction, delay differential equations. 



1. Introduction 

Activation cascades are frequently found in biological signal transduction 



systems [12|, |l8j . Perhaps one of the best studied examples is the mitogen- 



activated protein kinase (MAPK) cascade, which plays a central role in im- 
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Figure 1: (Colour) A cascade of length n. The nodes m the cascade can either be in an 
inactive state a;,*, or active Xi, shown graphically in the nodes of the right by the addition 
of a phosphate group. An external signal R{t) activates the first node. Once a node is 
active, it can activate the next node, and so on until the end. The activation rates are aj 
and the inactivation rates are Pi. Image adapted from 12 1. 



portant cellular functions such as regulation of the cell cycle, stress responses 
and apoptosis [isj . In general, activation cascades are formed by a set of com- 
ponents (typically proteins) that become sequentially active in response to an 
external signal (see Fig. [T]) . The role of cascades is to relay, amplify, dampen 
or modulate signals in order to achieve a variety of cellular responses. 

Activation cascades, particularly the MAPK cascade, have been the sub- 
ject of numerous studies, experimental and theoretical jo], 0, Isl, [l2[ 15, 16 



23l . |26| . In this work, we study ODE models of linear cascades and obtain 
analytical solutions in terms of the lower incomplete gamma function for the 
case when inactivation rates are identical and study the case in which a single 
protein has a different inactivation rate to the rest. Finally, we discuss how 
these results may be used in parameter fitting and model reduction as an 
alternative to delay differential equations. 

2. Linear cascades and their gamma function representation 

Consider a cascade of length n subject to an external signal R{t). When 
receiving the input -R(t), the first inactive component [x]) is transformed 
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into its active state (xi) which then activates the next inactive component 
(X2). Sequential activation of x* by continues until the end of the cas- 
cade. The output of the cascade is the active form of the last protein, Xn- In 
the case of the MAPK cascade, the components are proteins and the activa- 
tion corresponds to a post-translational modification, i.e., phosphorylation. 
However, the formalism below could also be used to describe other sequential 
biochemical processes with similar functional relationships, e.g. n-step DNA 
unwinding 17 . 

The mass-action kinetics model that describes the time evolution of the 



activation cascade is 12 



dxi 

dX2 

'df 



R{t){Ti - xi) - f3ixi, 
a2Xi{T2 -X2)- P2X2, 



(1) 



dXn 
dt 



where (Tj — Xi) = x* is the inactive form of Xj when the total amount of each 
component is given by Tj = Xi + x*. We assume resting initial conditions 
(i.e., Xj(0) = 0, for all i) and that Tj remains a constant (i.e., no protein 
production). As shown in Ref. [I^, whenever Tj 3> Xi we have Ti — xi ^ Ti 
(the so-called weakly- activated case) and we can re- write the system ([1]) as 



dxi 

"dT 
dx2 

~dd 



Rit) - Pixi, 

a2Xi - (32X2, 



(2) 



dXn 
dt 



where = djTj and R{t) = R(t)Ti. 

The system of equations ([2]) is linear and can be written in vector form: 



X = Ax + R(t)ei, 



(3) 
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where x = [xi, 



, Xn]"^, the n X n rate matrix A is: 



A 



-Pi 
a2 



-(32 



a. 



-(3n 



(4) 



and ei = [1, 0, . . . , 0]^ is the first n x 1 vector of the canonical basis. In 
general, we use Cj to denote the i-th canonical vector. 

2.1. Constant stimulus 

In an experimental setting, one often wants to study the response of a 
biological system to a constant stimulus (e.g., constant temperature, light or 
treatment): 

R(t) = ai G M, t>0. 
Then the solution to equation ([3]) with initial condition x(0) = is: 

x(t) = aiA-i[e*^-I„]ei, (5) 

where I„ is the n x n identity matrix, and e*^ is the matrix exponential. 

Recently, it has been shown that, given a fixed gain, a linear cascade 
provides optimal amplification when all the inactivation rates are identical 
(i.e., Pi = (3 for all i) and the cascade has a finite length [3]. 

Assume that our cascade is optimal with /3j = /3, then the rate matrix 
becomes 

^ -P 
a2 -P 



ar. 



-P 



and it can be shown (see Appendix A.l and Ref. [17 
that the solution for the output of the cascade is: 

x„(t)= (^yP(n,/3t), 



(6) 



for detailed calculations) 



(7) 



where P(?7,,/3t) is the normalised lower incomplete gamma function and a(„) 
is the geometric mean of the activation rates 

l/n 



a 



in) 



n 



(8) 
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2.2. Exponentially decreasing stimulus 

When the first protein in the cascade is subject to an exponentially de- 
caying stimulus (e.g., when the signal is a reactive molecule or it becomes 
metabolised, or the receptors become desensitised) 

R(t) = aie-^\ 

then the solution to Eq. ^ with initial condition x(0) = is 

x(t) = ai [e*^ - e-^*I„] A'^ [l„ + AA-^] ' e^. (9) 

If we assume again that (3i = P for all i (i.e., A = A), then the output of the 
system is given by: 

f te)%-'*P(r^,(/3-A)t) if/3^A 

Xn{t)=l (10) 

[ r(iT) if /3 = A, 

where is defined in ([8]) and r(n + 1) is the gamma function. As for the 
case of constant stimulus, the solution is also given in terms of the lower 
incomplete gamma function (see Appendix A.l for calculations). 



3. Perturbation of a single inactivation rate 

We now examine how the output of a weakly-activated (linear) activation 
cascade is modified when a single protein in the cascade has a different in- 
activation rate. For instance, Ref. [7| considered a cascade with an auxiliary 
protein with different inactivation rate (i.e., a leaky integrator) at the end of 
the cascade. We study the effect of such a 'perturbation' and how the effect 
depends on the position of the component in the cascade. 

Consider a cascade of n proteins with activation rates aj and inactivation 
rates (3j = 13, Vj ^ i, and Pi = (3 + e for a given node i. In this case, the rate 
matrix of the system ([3]) is 

A = A-eSieJ = Hi, (11) 

where A is given in Eq. ([6]) and ej is the i-th vector of the canonical basis in 
M". 
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The Jordan decomposition of Hj is 



Hj = QjJQj \ 
where J is the Jordan normal form 

-P 1 



(12) 



1 

-/3 



(13) 



and Qj is the matrix with generahsed eigenvectors as columns. 

An important property of this Jordan decomposition is that both J and 
the vector Q~^ei are independent of the location of the perturbation i. As 
shown in Example [6l this follows from the following fact: consider i < h 
(without loss of generality), then rows 1 to i — 1 and h to n of Qj and are 
identical, i.e., Qi(l : i — 1,:) = Q/i(l : z — 1, :) and Qi{h : n, :) = Qh{h : n, :) 



in Matlab notation. (See proof in Appendix B 



3.1. Constant stimulus 

The constant stimulus solution ([5]) for the perturbed case (A = Hi) is 



x(t) = «iH-i [e*"' - 1„] ei = aiQ,J-i [e'' - I„] Qr'ei, 



(14) 



which follows from ( fT2|) . The properties of J and Qj stated above imply 
that the vector [e*"' — 1„] Q~^ei is independent of the position of the 
perturbation, i. Hence the entries of x(t) are determined only by the matrix 
Qj. However, for two cascades with modified decay rates at i and h {i < h), 
we have Qi{h : n, :) = Qh{h : n, :), meaning that the solution for the last 
n — h + 1 proteins is the same in both cascades. (Similarly, the i — 1 first 
components of the solution ( IT^ are identical, but that is obvious.) Hence 
the position of the perturbation in the inactivation rate of a component does 
not affect the final output of the cascade. 

Example 1. Consider two cascades of length n = 6 with constant stimulus 
and activation rates ai = ■ ■ ■ = ae = 3 and degradation rate /3 = 2 for all 
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Figure 2: (Colour) Time course solutions of two e-pcrturbed cascades from Example [T] 
A: Cascade with a perturbation in the degradation of the third protein. B: Cascade with 
a perturbation in the fifth protein. The activity of the proteins of both cascades in the 
same for nodes above and below the perturbations (continuous lines), but is different in 
the nodes between the perturbations (dashed lines). 



proteins except for a perturbation e = 0.5 on the third and fifth proteins of 
each cascade, respectively. The corresponding rate matrices are: 



-2 

3 -2 

3 -2.5 

3 -2 

3 -2 

3 -2 

The Jordan form for both cascades is: 

-2.5 



-2 

3 -2 

3 -2 

3 -2 



-2 1 



-2 1 

-2 1 



-2 1 



3 -2.5 
3 



-2 
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and the corresponding generalised eigenvector matrices are: 



Q3 = 



Q5 = 
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(15) 



As explained above, rows 1-2 and 5-6 of Q3 and Q5 (in bold) are the same. 
In addition, 



Q3 ei = Q5 ei = 



Since the rows of of Q3 and Q5 are the same below the second perturbation, 
then the values of x^it) and XQ{t) are equal in both cascades. Figure |2] shows 
the time course of the proteins in both cascades: Xi(t), X2(t), x^it), and XQ{t) 
(solid lines) are the same in both cascades, while xsit) and X4^(t), the proteins 
"sandwiched" between the perturbations (dashed lines), are not. We discuss 
an application of this result in Section 14.21 
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1 
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Figure 3: (Colour) Time eourse solutions of two e-perturbcd easeades with decaying stim- 
ulus from Example [2j A: Caseade with a perturbation in the third protein. B: Cascade 
with a perturbation in the fifth protein. On the right of each plot, the label of each protein 
is placed at the level where its solution peaks. The activity of the proteins is the same 
for nodes above and below the perturbations (continuous lines) , but different in the nodes 
between the perturbations (dashed lines). 



3.2. Exponentially decreasing stimulus 

Just as in the previous section, the solution (|9]) for the exponential stim- 
ulus in the perturbed case (A = Hi) can be rewritten as: 

x(t) = «iQ, [e'' - e-^X,] J ^ [Irr + XJ-'] Qr'ei, (17) 

and, again, the same argument follows to conclude that the last n — h + 1 
components of the solution ( 1171) are the same for two cascades modified at 
positions i and h {h> i). 

Example 2. Consider the same cascades as in Example [1] with stimulus 
R{t) = aie~^* with A = 1. Figure |3] shows the time behaviour of the proteins 
in the two cascades. As in the previous example, the proteins above and 
below the perturbations are unchanged. 

4. Applications 

4.1. Model simplification and parameter fitting 

The expressions for the output of the cascade x„(t) in terms of incomplete 
gamma functions can be useful to fit activation data to a reduced number 
of parameters. Rather than fitting the observed output of the cascade to 
n + 2 parameters for an entire module with n components, the approximate 
expression with the lower incomplete gamma function contains at most four 
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Figure 4: (Colour) A: Schematic of a signal transduction model with a cascade. Given a 
linear cascade of arbitrary length, the red node node at the top is the stimulus, and green 
nodes are the components of the cascade. The last node of the cascade transmits the signal 
to downstream components of the pathway. The model of the cascade has up to n + 2 
parameters: ai, . . . , an, /?, and if the stimulus decays, A. Then the cascade is condensed 
into an expression with an incomplete gamma function that sends the exact same signal 
as the cascade in the left panel directly to the rest of the network. The new expression has 
parameters Q!(„), /3, n, and if the stimulus decays, A. B: Examples of the time-course of two 
cascades with constant (top) and exponentially decaying stimulus (bottom). The dashed 
lines indicate the solutions to the corresponding systems of ODEs, squares are noisy data 
generated from the models, and bold lines arc fits to the data using the incomplete gamma 
function expressions (see Example [3]) . 
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parameters: /3, n (and in the case of exponential decay, A). In this 

approach, (shown graphically in Fig. SJ/^) the length of the cascade n is turned 
into a fitting parameter, similarly to what is done with Hill coefficients. 
Indeed, the fitted value does not need to be an integer because the lower 
incomplete gamma function P(n, t) is defined for any positive real number in 
its first argument [H. 



Example 3. Consider two cascades of length n = 5 with parameters ai = 3, 
ttj = 4 for i = 2, . . . , 5 (so = 3.776), and /3 = 3. One cascade is subject 
to a constant stimulus R{t) = ai and the other to an exponentially decaying 
input R{t) = aie~^^ with A = 1. 

After solving numerically the ODE models of the linear cascade (|3]) for 
both inputs (dashed lines in Fig. HP), we sample the output x^it) at times 
t = {0, 1, ... , 10} and we add random noise from a distribution J\f{0, 0.05^) 
to generate our 'observed data' (squares in Fig. HP). 

We fit the gamma function expressions ([7]) and ffTOj) to the 'data' using the 
parameter fitting method introduced in Ref. [ij (see Appendix D). The bold 



lines in Fig. |1J3 show the fits to both cascade output data. The fits to the 
noisy data are good and the estimated values are close to the 'true' ones: in 
the case of the constant stimulus cascade, the fitted values are ~ 4.068, 
(3 3.281, and n ~ 5.418; in the case of the exponentially decaying stimulus, 
the estimated values are «(„) ~ 3.317, f3 ^ 2.177, n ~ 4.6, and A ~ 2.177. 

4-2. Cascade equation reordering 

The results presented in Sec. |2] and |3] allow us to reshuffle equations of 
cascade models where perturbations are known to occur. In particular, all 
proteins with the same inactivation rates can be grouped together upstream 
in the cascade so that they can be replaced with the incomplete gamma 
function expression, while the perturbed proteins are placed downstream 
and take the gamma function as an input. 

This process of reordering the cascade, which is schematically represented 
in Fig. OA., can be used to reduce the ODE model for the cascade without 
altering the dynamics or the timescales. Suppose we have an e-perturbed 
cascade of n + 1 proteins that we have reordered so that the first n proteins 
have inactivation rate (3 and the (n + l)-th protein has rate (3 + e. 
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Figure 5: (Colour) Example of cascade reordering and substitution. A: A linear 
e-perturbed cascade model of length 4, the input (red node) can either be constant or 
decaying. Green circle nodes are proteins whose inactivation rates are all /?, the blue star 
node has inactivation rate (5 + e. Downstream of the cascade lie other components of the 
signalling pathway. Reordering of the equations: the protein with the perturbed inactiva- 
tion has been moved the bottom of the cascade. Both this cascade and the one on the left 
have the same output. The first three equations in the reordered cascade are substituted 
for an incomplete gamma function. B: Numerical example of cascade reordering. Top: the 
equation of the perturbed protein is placed at the bottom of the cascade (the time-course 
of the untouched cascade is shown in Fig. . Bottom: the solution to the module of 
unperturbed proteins is given by equation (jlOp ; the solution of the perturbed protein at 
the bottom of the cascade is given by equation ((2T|) (see Example |4|). 



4-2.1. Constant input 

Consider first a constant input R{t) = ai. Using 
namics of tlie output of tlie cascade as 

dXn+l 



dt 



a 



/3 



, we write the dy- 
P(n,/3t)-(/3 + £)x„+i. (18) 



This equation can be solved analytically (see Appendix C.l): 



an+1 ( "(n) 



f3 + e V f3 



(£""^--(-^)"-*^)(/3t) 



k=0 



(19) 
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where we have taken the initial condition x„_|_i(0) = 0. 



4-2.2. Exponentially decaying input 

Consider an exponentially decaying input R{t) = aie~^^. First assume 
(3 X, and use ffTOl) to write the ODE for the {n + l)-th protein: 



dXn 



+1 



dt 



a(n) 



-xt 



P{n,{(3-X)t)-{f3 + e)xn+i. (20) 



When the initial condition is x„+i(0) = 0, the analytical solution for this 
equation is (see Appendix C.2 ): 



a I, 



(n) 



/3- X + e \f3- X 

n 

fc=0 



-xt 



—e 



(A - /3)"-^') (/3 - X)H'' 



(21) 



When /3 = A we have: 

dxn+i _ 
dt 



i n + 1 



(22) 



and the solution is 

Xn+l{t) - ' 



E 

A;=0 



■Iff 



k+n—k 



+ (-1 



(23) 



Example 4. Consider the n = Q cascade from Example Ej where the degra- 
dation rate of the third protein is e-perturbed (time-course shown in Fig.[2jA.). 
The equations of the system can be reordered without affecting its final out- 
put so that the equation of the perturbed protein is at the bottom of the 
cascade (Fig. |5j3, top). The output of the first 5 equations in the reordered 
cascade is then given by the incomplete gamma function expression fITU]) and 
the analytical solution of the perturbed protein (which is now the output of 
the cascade) is given by Eq. f l?T]) (Fig. [Sf3, bottom). 
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Figure 6: (Colour) A: Example of the use of linear activation cascades to replace delay 
differential equations. A signal node (red node) activates a node in a signalling pathway. 
The bottom node responds with a delay r. The delay in the equation is removed and 
substituted with a linear cascade of length n. The entire cascade is replaced by a lower 
incomplete gamma function. B: Top: The dashed line is the solution to the DDE (|24p from 
Example [Sj the squares are points taken from the solution with added random noise. The 
continuous line is the approximation using a lower incomplete gamma function. Bottom: 
the relationship between the ratio n//? and r is linear (dashed-line) with slope 0.977 and 
intercept 0.962. Inset: a almost does not vary with r (see text). 



4-3. Delay differential equation models for activation cascades 

Experimental observations in signalling cascades are typically concerned 
with the effect of the cascade on the output characterised in terms of the 
amplification, distortion and delay introduced in the output downstream. 
Within the framework of ODE models, the interactions between the vari- 
ables occur instantaneously. Hence, if the response to a stimulus occurs 
with delay, one must incorporate further intermediate variables that were 
not considered in the original ODE model, corresponding to unmeasured, 



hidden processes that take time to complete [21|. This process can lead to 



large models with many unobservable variables and large numbers of param- 
eters or to the introduction of abstract variables to model unknown processes 
that may contribute to the observed delays 0, 13|. Alternatively, modellers 
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often use delay differential equations (DDEs) to account for the wait between 
an event and its effect parsimoniously jsi ^ 19|. In a DDE, the activity of a 
variable depends on the state of the system in the past: 

| = fW«-r)). 

where the parameter r > is the delay. Linear systems of delay differential 
equations can be solved analytically using infinite series involving the Lam- 
bert function 0, 25|, but such solutions are often impractical to use. Our 



results indicate that simple delays can be modelled with linear activation cas- 
cades leading to concise ODE descriptions in terms of the lower incomplete 
gamma function and without relying on DDEs, as shown in Figure |6j\. 

Example 5. Consider a system with a delay, which we model with the fol- 
lowing linear DDE: 

dpi 

^ = dpi(t-r)-/3p2. (24) 

Figure [6|3 (top plot) shows the simulated time course of P2(^) (red dashed 
line) when ci = 2, /3 = 3, and r = 2 with initial conditions pi(0) = ^2(0) = 0. 
To generate our 'observed data,' we sample p2 at various time points and 
add observational random noise from a distribution A/'(0, 0.05^). 

We can fit this noisy data to a linear cascade of length n under constant 
input with parameters and (3: 

Pn{t)={^y nn.pt) ^p,{t), (25) 

and estimate the corresponding parameters. The solid line in the top plot of 
Fig. IHB shows the best fit of the data to a linear cascade, as obtained with 



a parameter fitting algorithm (see Appendix D). The estimated parameter 



values are a(„) ^ 2.27, (5 7.53, and n 22.1072. 

We also explore the connection between the parameters of the DDE and 
the best fitted (linear) activation cascade model, in particular as a function 
of the delay r. We simulate the DDE flMj) with parameters a = 2 and /3 = 3) 
for different values of the delay r G [0,5] and collect data from these models 
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as above, but without adding random noise. The dependence of the fitted 
parameters and r is shown in Fig.EP (bottom plot): a remains relatively con- 
stant while the ratio nj /3 grows almost linearly with t n/ = 0.962 + 0.977r. 
In fact, /3 increases linearly (slope ~ 4.11) and n grows almost quadratically 
with r (exponent ^ 1.88). When the delay r in the 'data' is increased, the 
length of the fitted cascade (n) increases while the time scale of the gamma 
function (1//3) decreases in consortium. If one attempts to fit the data al- 
lowing only the length fitting parameter, the fit is not successful, 
thus underscoring the importance of the time scale /5 in the approximation. 
Indeed, the original time delay r in the DDE is approximated in the lin- 
ear cascade by (n//3 — 1), i.e., the accumulated time needed to traverse n 
sequential steps with duration 1//3. 



5. Conclusions 

We have examined the classic model of activation cascades in the weakly 
activated limit 12 1. Following recent results [3], we have considered the 



important case where all inactivation rates of the components of the cascade 
are identical, which has been shown to provide optimal amplification. Our 
results show that the output of these cascades can be represented exactly by 
lower incomplete gamma functions. We also show that when one protein in 
the cascade has a different inactivation rate than the rest, its position in the 
cascade does not affect the final output. We have used these results to reduce 
the number of equations and parameters in ODE models without affecting the 
dynamics or the timescales of the system. We also show that in some cases 
incomplete gamma functions can be used to approximate delay differential 
equations. Beyond its application to enzymatic activation cascades, similar 
mathematical models of cascades could be helpful for the parametrisation and 
modelling of multi-step transcriptional processes, an area of active research 
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Appendix A. Calculation of the output of the linear cascade, Xn{t) 

Appendix A.l. The incomplete gamma function: Definitions and notation f^] 
The gamma function is defined as: 



oo 

a-1 



T{a) = / e-'s^-^ds, Re(a) > 0. (A.l) 
Jo 

One property of the gamma function is r(?7,) = (n — 1)!, n G N. 
The lower incomplete gamma function is given by: 

7(a,t)= [ e-'s"-Ms, Re(a) > (A.2) 



A different way of writing 7(0, t) when a = n G N is 20 



/ ,k\ 
7(n,t) = (n-l)! l-e"*5^M. (A.3) 
\ k=o ' J 

The normalised lower incomplete gamma function, which we use in our 
calculations, is defined as: 

P(a.O = ^. (A.4) 

r(a) 

Appendix A.2. Constant stimulus: derivation of Eq. ^ 

The solution of the ODE system ([3]) for x(0) = is given in Eq. ([5]): 

x(t) = aiA-i[e*^-I„]ei. (A.5) 

It can be shown by mathematical induction that the n-th component of the 
solution is: 



A similar result is obtained in Ref. [17| from the analysis of linear models of 
DNA unwinding. 

Using (1A.3P and the properties of the gamma function equation (1A.6P 
becomes: 

-nit) = (^y%^= ( ^Vp(n,/3t). (A.7) 



/3 ; Tin) V /3 
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Appendix A. 3. Exponentially decaying stimulus: derivation of Eq. IfT^) 
The solution of the ODE system ([3]) for x(0) = is given in iQ: 

x(t) = ai [e*^ - e-^X] A'^ [l„ + AA-^] ' e,. (A.8) 
When /3 7^ A, we can use mathematical induction to show that 



Xn{t) 



a 



in) 



(3-\ 

a(n) 
f3-\ 



n-l 



E 



t''{(3-\y 



k=0 

-xMn, (/3-A)t) 



(A.9) 



Tin) 



a 



(n) 



(3-\ 



e-^T(n,(/3- A)t). 



(A.IO) 



When /3 = A, we solve sequentially the ODEs ([3]) and use mathematical 
induction to get 



Xn{t) 



-/3t 



r(n + i) 



(A.ll) 



Appendix B. Properties of the e-perturbed matrix of rates, Hj 

The matrix Hj, which corresponds to a linear cascade with a perturbation 
e at position i, is defined in Eq. (1111) and has a Jordan decomposition given 
in Eq. ([12]): 

H, = QiJQr'. (B.l) 

As can be seen from the lower-triangular structure of Hj, its Jordan 
form is the direct sum of two Jordan blocks associated with —{(3 + e), with 
multiplicity 1, and — /3, with multiplicity n — 1: 



3 = 3. 



where 



J-{/3+e) = [-(/3 + £)]lxl' 



-/3 1 



(B.2) 



1 

-/3 



{n-l)x(n-l) 

(B.3) 

The matrix Qj contains the generalised eigenvectors {q^}"=i as columns |lO| : 

Q.= [qi|q^2l---|qy • 
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Proposition 1. The following properties hold for Qj and J: 
I. J is independent of i. 

II. q^, the eigenvector of Hj associated with —{f3 + e), is given by: 



qi(j) 







if j < i, 



(B.4) 



III. {q^n-k+i}k=i^ ('^ ~ -'-) generalised eigenvectors of Hj associated with 
—(3, are given by: 



If i < A; 

qn-fe+l(j) 

lii> k 

Qn-fc+lO) = 



(^) 



l\j-k 



a 



(k) 



i e {k + 1 . . .n}. 

j = k, 



(B.5) 



(B.6) 



IV. The inverse of the transition matrix Q ■ ^ has the following structure: 
• The first row is 



if j < h 



(B.7) 



if j > i- 



Rows e {2, . . . , n} are 

' iij<n-k, 



-~— ifj = n — k + l, 



ii j = n — k + 2 and n — i + 1 > k, 
otherwise. 



(B.8) 
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V. For all i: 



Qi = ei 



1 




(B.9) 



Proof. I. All the matrices Hj have the same eigenvalues with the same 
multiplicity so the Jordan form J is identical for all i. 
11. Let v} = Hjq]^, if i = 1 it is straightforward to see that v];(l) = — 
and when j > i 



i-2 



' f3a' 



(i) 



so ^rl = -{/3 + e)q\. 

When i > 1 then "v\{j) = for j = 1 . . . i — 1. U j = i, then 



(B.IO) 



vi(z) = -{(3 + e)a\i) [ 



i-l 



and if J > z then the same situation as in equation fIB.lOp applies, so 
again we have = — (/3 + £:)q^. 
III. Define Bj = (Hj + we can see from the definition of q2 in equa- 
tions fIB.Sp and flB.6p that BjQg = 0, i.e., is the eigenvector of Hj 
associated to —f3. The rest of the generalised eigenvectors qg . . . qj^ 
associated with — /3 can be multiplied by Bj to show that 



h = 3 . . .n. 



(B.ll) 



It follows that Bf = 0, so Qj as defined in equations (lB.4p . ( IB.Sp . 
and (]B.6P is the transition matrix of generalised eigenvectors of Hj. 
IV. We can verify that defined in equations in equations (1B.7P and 
(IB.Sp is the inverse of the transition matrix Qj by multiplying them to 
obtain Qj~"'^Qj = In- 
V. This property follows directly from the structure of Q~^. 



□ 
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Example 6. Consider a cascade of length n = 6. 
transition matrix and its inverse for i = 1, . . . ,6 are 



The structures of the 
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^3 
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(3) 
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(6) 




2 

^3 
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,5 
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^(1) 

^2 
(2) 



1(1) 



,3 
(3) 



,5 
(5) 



1(1) 





^4 
(4) 



,5 
(5) 



Qi 



'(1) 





"(1) 








^(1) 








^(1) 





q4 = 



(2) 



(2) 




(2) 




(2) 




(3) 




(3) 




(3) 




(3) 




(3) 




(3) 




(4) 




(4) 




(5) 










(B.12) 
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-f— 

'(4) 
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(B.14) 
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(B.15) 
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(B.16) 



(S) 








(B.17) 



Appendix C. Calculation of Xn-\-i{t) with one e-perturbed inacti- 
vation rate 

Appendix C.l. Constant stimulus, derivation of Eq. / fTPj) 

We solve the differential equation f|T8|) through an integrating factor to 

get: 

e(^+^)*x„+i(t) = {^y-J I e('^+^)*P(n,/3t)dt + c. (C.l) 
Use the properties of the Gamma function to re-express the integral as 

f e(^+^)*P(n, /3t)dt = f e(^+^)*dt - f e''Yl ^^^^ (^.2) 

and solve the second integral of (1C.2I) using integration by parts: 

fe=0 fc=0 \fc=0 ^ ' / 

(C.3) 
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to obtain 



1 + ^ 

e 



n-l 



A:=0 



k\ 



+ 



e — ' k\ sin — 1) 



The integral on the right-hand side of equation ( 1C.4P can be solved using the 



formula 1 1 



.fc=0 



e^+\n -\-k)\ 



n—\ — k 



(C.5) 



Substituting in equation flC.2p and gathering terms we obtain 



e(^+")*P(n,/3t)dt 



whence Eq. (IC.ip becomes 



k\~^ 5^+i(n- 1 - A;)! 



fc=0 



+ c, 
(C.6) 



(3 + e V /3 



n-l 



fc=0 



k\ efc+i(n- 1 - A;)! 



-|-ce' 



The initial condition a;„_|_i(0) = requires that 

^ ' f3 + e \ e 
which gives the final expression given in Eq. ([19 



(C.7) 
(C.8) 



Xn+l(t) 



/3 + e \ (3 



1 -e" 



k=0 



Tn-kJ^\ 



(C.9) 



Appendix C.2. Exponentially decaying stimulus, derivation of Eqs. I[21\) and (d^j 

Consider first the case P ^ X. To solve the differential equation fl2U]) . 
define a = /3 — \ and use integrating factors to get: 



(v')"(/~-/^-"|:^*)+-«=.io) 
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Using integration by parts for the second integral on the right-hand side, and 
following similar steps to those above gives 



a + e \ a 



fc=0 



The initial condition x„+i(0) = requires that 



a + e \ e 



+ c. 
(C.ll) 

(C.12) 



thus giving equation (12T|) : 



n-1 



g-(/3+e)t 



E 

fc=0 



(C.13) 



When A = /3 we use integrating factors to solve equation (l22il : 



a;„+i(t)e 



r(n + i) 



re«dt + c 



r(n + 1) e ^ e'=(n- A;)! 



+ c, 



(C.14) 



and taking the initial condition a;„+i(0) = we get Eq. f l23|) : 



n 



\kj-n~k 



A;=0 



Appendix D. Parameter fitting 

We use an evolutionary Monte-Carlo optimisation method with local 
search acceleration to fit the parameters of our models to 'observed data' 
in Examples |3] and [5l The error function to be minimised is a quadratic 
distance (sum of squared errors) between the data and the fitted values. In 
all the fits, the optimisation algorithm was run starting from an initial prior 
distribution f/(0, 10) for each parameter and the algorithm was run for 10 it- 
erations sampling 500 points at each iteration and keeping the best 50 points 
as a prior for the next iteration. The method is described in detail in Ref. jsj. 
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